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Abstract 

We continue the investigation on the applications of the Kramers equation 
to the numerical simulation of field theoretic models. In a previous paper we 
have described the theory and proposed various algorithms. Here, we com- 
pare the simplest of them with the Hybrid Monte Carlo algorithm studying 
the two-dimensional lattice Gross-Neveu model. We used a Symanzik im- 
proved action with dynamical Wilson fermions. Both the algorithms allow 
for the determination of the critical mass. Their performances in the defi- 
nite phase simulations are comparable with the Hybrid Monte Carlo. For the 
two methods, the numerical values of the measured quantities agree within 
the errors and are compatible with the theoretical predictions; moreover, the 
Kramers algorithm is safer from the point of view of the numerical precision. 
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I. INTRODUCTION 



In a previous work [|I|, we used a compact operatorial formalism to study the appli- 
cability of the Kramers equation to the numerical simulation of quantum field theories on 
the lattice. In the standard Monte Carlo approach, we must generate ensembles of states 
distributed according to a given statistical weight. In our strategy has been to realize a 
diffusive process involving the physical degrees of freedom but also auxiliary variables which 
drive the diffusion. The final equilibrium distribution approximates the desired one. The 
resulting schemes can be made exact by a global accept-reject test and provide thus gener- 
alizations of the usual Hybrid Monte Carlo and of the exact hyperbolic algorithm P|. 
For these algorithms, a good behaviour with increasing volume and from the point of view 
of numerical precision is expected. Moreover, in [Q] we have shown how the freedom in the 
auxiliary variables sector may be used to reduce significantly the autocorrelations in the 
free field case. In this paper we study the behaviour of the simplest among the proposed 
schemes. We take the Hybrid Monte Carlo as a reference point. Our theoretical labora- 
tory is the two-dimensional A^-fermion Gross- Neveu model with Wilson fermions This 
model is an asymptotically free purely fermionic model with a discrete 75 symmetry whose 
dynamical breakdown leads to mass generation. It has an exact large solution and is 
expandable ^ allowing for analytical investigation of its non perturbative properties. 
In our lattice simulation of the model we use the Wilson formulation to avoid the doubling 
problem and let the fermion number N be free. However, using Wilson fermions, a bare 
mass must be introduced and tuned in order to restore the chiral symmetry in the con- 
tinuum limit. The restoration of chiral symmetry is a significant test for new algorithms: 
previous works showed that non exact algorithms (Langevin, pseudofermions) fail to 
locate the critical point by means of mixed-phase techniques. In this work we shall show 
that the algorithm based on the Kramers equation does not encounter this problem. The 
use of two-dimensional four-fermions models as prototypes for realistic theories like QCD is 
somewhat limited. However, a non trivial model with a discrete chiral symmetry breaking 
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and many dynamical fermionic degrees of freedom implies enough problems to be a severe 
test for our algorithms. 

In Sec. |T| we introduce the Gross-Neveu model. In Sec. |T| we describe explicitly the 
algorithms that we have used. In Sec. ^ we discuss how we compared them. In Sec. [V| 
we review the previous simulation of the Kramers algorithm and of the Symanzik improved 
Gross-Neveu model. We then describe our explicit numerical simulation and present the 



results. Finally, Sec. ^ is devoted to a short summary and conclusions. 



II. THE GROSS-NEVEU MODEL 

The action of the model in the continuum is 

^ = yrf\(^^(")^V^(")-i/(V^(")V^("))%mV^("V^"^) a = l---N. (2.1) 

where ip^"'^ is a multiplet of two-dimensional Dirac fermions. 

Rewriting the quartic interaction with a Lagrange multiplier we obtain 

S = j (fx + (f^ + m)^("V^"^ + ^^'^ • (2.2) 
The manifest U{N) symmetry of the model can be enlarged to 0{2N) by writing the 



fermionic fields in terms of their hermitean Majorana components. We refer to ||TOl for 
a discussion of the renormalization properties of the continuum model in the two formula- 
tions with and without the auxiliary a field. In the chiral limit the model enjoys the discrete 
symmetry 

^ 75^ -^75 cr -a. (2.3) 

As a consequence, the potential of a , related to that of the composite field ipip, is symmetric 
under the exchange a —a and possesses two degenerate minima. The chiral symmetry 
is spontaneously broken. The same holds on the lattice, but in the Wilson formulation we 
need an explicit mass term to control mass renormalization. 
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The expansion of the model is obtained by integrating out the fermion fields. By 
introducing the large- coupling A = Ng"^ and the field S = a + m, the corresponding 
"effective" action is 

S^ff = N j d^x(^^{T.-mf -Ti\ogM{i:)^ M{T)=^+T.. (2.4) 

In a perturbative expansion in powers of g, the fermions remain massless in the chiral 
model. However, already in the leading order of the expansion, a non-zero expectation 
value of the S field is dynamically generated playing the role of a non-perturbative fermion 
mass. From the study of the phase diagram of the model we can infer the existence of 
a critical point in the weak coupling region. At this point the chiral symmetry can be 
recovered (and dynamically broken) with an appropriate choice of the bare critical mass. 
The vacuum expectation value (S) is exponentially depressed as A — following the leading 
renormalization group prediction (S) = ±Aexp(— vr/A). Actually, on a finite lattice at non 
zero A one chooses the perturbative mass in order to have a double well shaped potential 
for the S field having two degenerate minima. The E — S symmetry is recovered only in 
the chiral limit; this is reminiscent of the doubling problem since one has exactly 

(F(7/>,V',S,r)) = (F(75^,-^75,-S,-r)), (2.5) 

where r is the Wilson (or Symanzik) coefficient of the mass term which solves the doubling. 
Quantities like the asimmetry in the steepness of the potential in the two minima can be 
1/A^ expanded and computed. 

Since we are in two dimensions, finite volume effects are particularly dangerous. There- 
fore we improved the lattice Wilson action at the tree-level following the method of 
Symanzik [^]. The resulting action improved at order O(a^) is 

- r 2 1 



T 



+ (2.6) 



E 



'ipni'nim + (7„) + ^cr^ 



where we have understood the flavour indices. The Fermi fields have antisymmetric bound- 
ary conditions in the temporal direction. The parameter is the Wilson parameter for the 
Symanzik action. We have taken = 1/3 remarking that the strong coupling matching 
between Wilson and Symanzik actions is at = 1/3 r^^. In the interested reader can find 
a detailed account of the leading and next-to-leading 1/A^ expansion of the model described 
by Eq. (p.6|) . In momentum representation, the fermionic matrix is 

M = i^^^p^ + Ms{p) + E, (2.7) 

where we have defined 

p^ = 2sin^, (2.8) 

On a finite L x T lattice, the form Eq. ( |2.4| ) of the effective action together with Eq. ( p.7| ) 
allows for the computation of any relevant quantity in the leading 1/A^ limit. 



III. DESCRIPTION OF THE ALGORITHMS 

We now describe in an unified scheme the algorithms that we have used. We give all 
the details showing how the fermionic fields are dealed with. Because of the particular 
structure of M(S), we may introduce real auxiliary fields x''"^ (ce = 1. . .N and Dirac 
indices understood) and conjugate bosonic momenta it to form the extended action 



V V 2A 2 2-^ 

sites \ 



M{i:fM{T) 



X 



(a) 



(3.1) 



For each given configuration of S and of the auxiliary pseudofermionic fields we define 
the following fields 



$(")(S) = [m(S)^M(S) 



X 



(a) 



(3.2) 
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where the sparse matrix inversion is performed by means of the Conjugate Gradient algo- 
rithm 

The X fields action is quadratic. Therefore, we can extract randomly the x fields accord- 
ing to their exact equilibrium distribution and evolve the other fields between two successive 
such refreshments. Note that this procedure introduces another free parameter, namely the 
frequency of the x fields update. 

The single sweep updating is given by the following algorithm 

• Refresh the auxiliary x*^"^ according to 

y(°) = M^(E) 

(3.3) 

TT = 7]. 

Here, t]^"-* and rj are vectors of gaussian random numbers with zero mean and unit 
variance. This step must be done only every k sweeps. 

Update the momenta vr according to 



r = e '"TT + Vl - ^, (3.4) 



r' — ^-^^n 



where ^ is a gaussian random number with zero mean and unit variance. The param- 
eters 7, e are positive real numbers. 

Integrate the equations of motion associated with the force 

F„(S) = -^(S„ - m) + i $(°)-^ (m^(S)M(S)) (3.5) 

{n is a site index) using Nmd iterations of the leap-frog scheme 

7r(t + e/2) = 7r(t) +e/2 F(S(t)), (3.6) 
S(t + e) = S(t) +£7r(t + e/2), (3.7) 
7r(t + e) = Ti{t + e/2) + e/2 F{T.{t + e)), (3.8) 

or its dual form in which the roles of S and vr are interchanged. 



• Perform a Metropolis test between the states before and after the integration of the 
equations of motion. If the test fails, reject the proposal for the S field and negate all 
the momenta 

S Soid TT -'"'old (3.9) 

otherwise accept the proposal for both S and vr. 

We have four free parameters: k, Nmd, e and 7. Detailed Balance is exactly satisfied for 
each set of values they take. The algorithm under study (which we shall call "Kramers 
algorithm") is obtained with N^d = 1 and arbitrary values of k, 7 and e. The usual Hybrid 
Monte Carlo corresponds to = 1, 7 ^ 00 and arbitrary values of N^d, ^- Of course, in 
the Hybrid Monte Carlo case, there is no mixing between the old and new momenta and it 
is unnecessary to negate them on rejection. 

We remark that from a physical point of view, the Kramers algorithm possesses two 
distinct time steps, ei„ = £7 drives the "irreversible" motion, whereas erev = ^ drives the 
"reversible" one. Apart from effects coming from the negation of momenta, we expect that 
the Kramers and the Hybrid Monte Carlo algorithms to behave similarly when e^^r ~ '^/Nmd 

with erev = ^md- 

IV. COMPARISON 

In |]T]] we argued that the Kramers algorithm should behave better than the Hybrid Monte 
Carlo for large volumes. The argument is rather naive and must be confirmed numerically. 
Consider a lattice model satisfying L ^ ^, where L is the lattice size and ^ the correlation 
length. On general grounds, the acceptance rate of the Hybrid Monte Carlo is known to 
be i 

Pace ~ eifcicNmde^Vv). (4.1) 

Optimal tuning is expected to require N^d ~ 1/e giving a sweep-sweep correlation r ~ 1/e. If 
the volume is varied, we have to scale e ~ V^^^'^ in order to keep the acceptance probability 
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constant. The advantage of the Kramers algorithm is that N^d = 1- If we neglect the 
influence of the momenta negation on the acceptance, then the scaling relation is modified 
to e ~ 1/-^/^ If, moreover, the optimal autocorrelations stay ~ 1/e then it follows that on 
large volumes the Hybrid Monte Carlo behaves worst than Kramers. 

Empirically the maximum acceptable value for e in the Kramers algorithm was rather 
greater than the one of the Hybrid Monte Carlo. This is reasonable because, at fixed e, as 
Nmd is increased the acceptance rate decreases before reaching a plateaux. 

Another advantage of the Kramers algorithm is related to the numerical precision, the 
opportunity of having Nmd = 1 without penalty is surely welcome; for instance, in QCD, 
as the volume (L^) grows or the quark mass (m) is decreased we must have N^d ~ L/m^^^ 
and some protection seems necessary to protect irreversibility against the accumulation of 
numerical errors. 

To compare quantitatively the performances of the two algorithms we must determine 
which is the computer time needed to have a given statistical error. If we denote with tq the 
integrated autocorrelation of the observable Q, then the computational cost of an algorithm 
measuring Q is proportional to 2tq^ which is the reduction factor which determines the 
actual number of independent measures. This factor must be multiplied by the number 
of Conjugate Gradient inversions needed to make a single sweep. Indeed, in a model with 
dynamical fermions it is reasonable to neglect the time expended during all the other steps 
of the algorithm. In our scheme we have blocks each made of N^d molecular dynamics steps 
followed by a Metropolis test and we refresh the pseudofermionic fields x only at the end of 
k such blocks. Let us call LFi and LF2 the two dual second order leap-frog schemes (see 
for a compact review of symplectic integrators) which we sketch as follows 

LFi : 7r(to) ^ S(t + e/2) ^ 7r(t + e), (4.2) 
LF2 : S(to) ^ vr(t + e/2) ^ S(t + e). 

Keeping x fixed, a new Conjugate Gradient inversion is needed only when the S field con- 
figuration or the X fields change. The advantage of the LFi scheme is that it computes the 
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old and new action at the same evolution times at which the force itself must be evaluated. 
The quantity determining the computational cost is the number of inversions per sweep 
computed between two refreshment of x 

C{k, Nmd)LF, = Nmd + 1/k, (4.3) 
C{k, Nmd)LF, =Nmd+l + l/k, 

and the actual performance is just this {Q independent) factor corrected by the [Q depen- 
dent) autocorrelation. From Eq.( |4.3| ) we see that LF2 must certainly be ruled out at small 
Nmd- For what concerns k, in the Hybrid Monte Carlo case the l/k term is usually not 
very important because N^d ^ 1- However, with the Kramers algorithm, a large k is twice 
better than k = 1. Moreover k > 1 may allow for a more efficient exploration of the phase 
space. 

We remark that our conclusions hold provided that LFi and LF2 have the same behaviour 
concerning autocorrelations. 

V. SIMULATION AND NUMERICAL RESULTS 
A. Review of previous results 

1. Previous study of the Kramers algorithm 

In , the Kramers algorithm was applied to the simulation of quenched compact QED 
on a 8^ lattice in the disordered phase. Measures of the average plaquette and of the average 
2x2 Wilson loop were made using also the Hybrid Monte Carlo algorithm as a reference 
method. The performances of the two algorithms turned out to be comparable. 

2. Previous study of the improved Gross-Neveu model 

In [Q, the Symanzik improved Gross-Neveu model was studied on the lattice using the 
pseudofermion algorithm. In this scheme the variation of the fermionic contribution to the 
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action is expanded to the lowest order in the S field variation ST,. 

5S~^ ^ {T- 5T + ^ 6i:^ -N TrM-^(E) ST. (5.1) 

The fermionic matrix inversion was performed computing 

M-,' = {x^x]) (5.2) 

as a Monte Carlo average in presence of the quadratic action = x^M~^x- The systematic 
error to be kept under control has two contributions. It gets an easily controllable statis- 
tical contribution from the inversion and a difficult contribution due to the truncation of 
the fermionic term at small but non zero 6T. An extrapolation as ((5S)^) must be 
performed. Using the mixed-phase tecnique, the critical bare mass rric was determined at 
different values of the parameters and of the step ((5S)^). The extrapolation did not agree 
with the theoretical prediction in the case of (E) and was explained as a wrong determi- 
nation of rric. Moreover, deviations from the Schwinger-Dyson equations of the model did 
not show the expected linear dependence on ((5S)^). These problems were absent for (S^) 
which is expected to have a milder dependence on rric. 

In [^], the same kind of measurements were done using two different algorithms: the 
Langevin one with bilinear noise and the Hybrid Monte Carlo. In both cases, the inversion 
of the fermionic matrix was performed by the Conjugate Gradient algorithm. The Langevin 
algorithm is not exact, given a step e the update of the bosonic field is computed as 

S S 

T{x, t + e) = E(x, t) - 6 + ri{x, t) + (5.3) 

dT[x, t) 

y,z,w OZ^{X,i) 

where rj and rif are gaussian fields normalized as 

{7]{x,t)r]{x',t')) = 26,,,6u'. (5.4) 

An extrapolation as e — is needed. The Langevin algorithm was found to be completely 
unable to determine the critical mass since the critical point was so unstable that different 
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seeds of the random number generator gave different results. Using the theoretical value for 
the critical mass the simulation in the definite phase was stable only for the smallest values 
of e. On the other hand, the Hybrid Monte Carlo simulations turned out to be perfectly 
stable in both the mixed and definite phase simulations. The results (included the stable 
Langevin runs) were consistent with the theoretical predictions. 



B. Our results 

In this work, we have repeated the measures of |^,^ using our favoured algorithm. More- 
over we have studied the variation of the performance as the tuning is varied, an information 
which is very important from a practical point of view. 

The numerical simulation has been done on the APE supercomputer ||15|. The model 



operating at Pisa is the so called "tube" machine, a 128 processor parallel computer with a 
peak performance of 6 GigaFlops. All the code has been written in the high-level "apese" 
language 



We have used the Conjugate Gradient algorithm (see [T^IT^ for a discussion oriented to 
the Dirac lattice operator) to obtain iteratively the inverse of the sparse symmetric matrix 
A = M'^{a)M{a). We did not use any kind of preconditioning, the improved Symanzik 
action has next-to-neighbours interactions and an efficient incomplete factorization is not 
trivial. Moreover, the condition number of A becomes larger when S ^ 0, but at A = 2.0 
on the 40^ lattice, the average S (which is the fermion mass) and its fluctuations provide a 
good rate of convergence and the inversion problem was not critical. We tried also to invert 
directly the matrix M using the Biconjugate Gradient algorithm [|l^ but we did not see any 
advantage. The stopping condition in the Conjugate Gradient algorithm is usually chosen 
to be a bound like 

\\r\\ < ecG\\b\\ (5.5) 

r is the residual vector: in the following we shall denote with r„ the recursive residual 
computed on fly by the algorithm and we shall denote with f„ the residual defined by 
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fn = 6 - AXn- (5.6) 

We stress an important fact: let x* be the exact solution, the Conjugate Gradient algorithm 
guarantees that \\xn — x*\\ monotonically decreases. On the other hand ||r„|| and ||f„|| 
may show relatively large fluctuations. Hence a decreasing trend must be checked on the 
succession of residuals to keep under control this effect. Let us now address the problem 
of numerical precision. In absence of rounding errors there would be no difference at all 
between the two residuals r„ and f„ as can be checked inductively using the explicit form of 
the algorithm. The effect of rounding errors associated to a finite precision machine can be 
easily studied on simple model problems (like harmonic problems on regular lattices). The 
result is that, given the exact solution x*, the recursive residual r„ tends to zero, while the 
true error has a positive infimum 



tij lY) tic 



\X^ 



e as n — > oo (5-7) 



of the order of the machine precision. The small number e may be further reduced if all 
the scalar products in the algorithm are implemented as binary sums. The same saturation 
phenomena occurs to the residual f, namely 

I |f„| I — >• 1 16| I as n — >• oo. (5.8) 

Here e' is a number which depends on the machine word length, but also on the norm of A. 
Typically the saturation of || *\\ and ||f„|| happens roughly at the same time. 

In a realistic problem x* is not known and one can study the behaviour of r„ and f„. The 
departure of Vn from f„ is typically very slow until the machine precision is reached and 
saturation sets up. Our stopping condition has been 

||r„||<10-^ (5.9) 

and we checked on configurations taken randomly during the runs that the saturation of 
I |f„| I was reached. 
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We used a 40^ lattice with = 10 flavours. In the framework of the 1/A^ expansion one 
can study the width of the scahng window, which is roughly measured by looking at the 
value of S in the crossover region where the weak and strong coupling expressions for the 
perturbative mass are equal. A detailed discussion of the scaling properties of the model 
including next-to- leading 1/A^ corrections and flnite size effects can be found in We have 
chosen A = 2.0 at which the correlation length is about three lattice spacings and we are in 
the scaling window. The volume factor in favor of the improved Symanzik action is about 
15 and at L = 40 we expect from the flnite lattice 1/A^ prediction to have very small flnite 
size effects. 

We determined the critical Wilson mass at which a flrst order transition takes place using 
mixed-phase runs and measured in the two deflnite phases (S) and the composite (S^)c which 
provides informations on the steepness of the vacua. For these quantities the next-to-leading 
1/A^ corrections can be computed and compared with the numerical simulation. 

On the 40^ lattice the theoretical predictions are rric = —0.974 and 

(S_) = -0.293, (S2), = 0.190, (5.10) 

= 0.311, (S^)c = 0.146. (5.11) 

We started taking k = 1. After a very rough tuning of 7 and e using deflnite phase runs 
with the 1/A^ theoretical mass, we choose 7 = 5.0 and e = 0.06 for the mixed-phase runs. 
Setting half of the lattice to the positive vacuum and the other half to the negative one, after 
a short thermalization, a kink conflguration appears which is very long-lived when the bare 
mass is at its critical value. From inspection of the leading effective potential one checks 
that bare masses smaller than the critical one favours the negative vacuum and viceversa. 
In Fig. (I) we show the time evolution of (S) for 16 different values of the bare mass starting 
from rric = —0.958 separated by Arric = 0.002. The result is 

= -0.969 ± 0.002, (5.12) 

In this part of the simulation we did not insist on a heavy tuning of 7, e and k because 
we were mainly interested in the exact determination of rric. We then turned to a detailed 
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tuning of 7 and e using definite phase runs at the critical point. During these runs k was set 
to the unoptimized value 6. In Tab. (|) we report measures of autocorrelations obtained in 
the two phases of the model with different values of 7 and the common choice e = 0.09. Pace 
was about 70-75% and greater values e did show a sharp decrease in Pace- We did not find 
any sistematic dependence of 7 for the two measured expectation values of (S) and (S^)^. 
All autocorrelations are computed with the automatic windowing algorithm by vary- 



ing the c constant in the range 4.0-8.0, cross-checking the result with the statistical ineffi- 



ciency [|^] and using consistently at least 1500r sweeps. 

In the mixed-phase, unoptimized Hybrid Monte Carlo runs gave the same value for the 
critical mass. We did not compare the performances because optimization in mixed-phase 
runs was too time consuming. In the definite phase runs we have taken e = 0.06 and 
Nmd = 8 corresponding to a ratio of C factors between the two algorithms of about 7.7. The 
related rescaled autocorrelations are shown in Tab. (|I|). The resulting expectation values 
were always consistent with the Kramers algorithm and the predictions. 

Best estimates for the average values extracted only from the runs with the Kramers 
algorithm at different 7 are 

(S_) = -0.279 ±0.002, (E2)^ = 0.189 ±0.001, (5.13) 

= 0.316 ±0.002, (S^)c = 0.145 ±0.001, (5.14) 

and are compatible with both the theoretical predictions and the Hybrid Monte Carlo results. 

From the data obtained we see that in the negative phase 7_ = 1.0 is the optimal 
value for both the operators. The same happens in the positive phase but with a smaller 
optimal value 7+ = 0.5. We can see that in the negative phase the performances of the two 
algorithms are comparable, whereas in the positive one the Kramers algorithm is about 1.5 
times slower. 

The existence of an optimal range for the parameter 7 at fixed e is easily understood as 
follows. In the 7 ^ 00 the Hybrid Monte Carlo with Nmd = 1 is recovered and it is well 
known that this is far from being an optimized regime. In the opposite limit, 7 — > 0, we 
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obtain the classical equations of motion in the background x fi^ld between two refreshments. 
Of course molecular dynamics without noise is highly self-correlated and again a large r 
results. In the Gross-Neveu model the 7-^0 limit is particularly dangerous. We can 
study the N = 00 model taking for rric the theoretical prediction including next-to-leading 
correction on the 40^ lattice at A = 2.0. One finds that the S effective potential 
is so unbalanced that one of the two minima disappears. In other words, if the quantum 
(noise) dynamics is switched off keeping all the parameters fixed, the behaviour of the system 
changes dramatically. 

VI. SUMMARY AND CONCLUSIONS 

This work is a first numerical application of a previous investigation on the Kramers 
equation approach to the lattice simulations of field theoretical models. We proposed many 
schemes among which the simplest (introduced in is studied here. The situation which 
we have chosen should be one of the worst possible cases, namely a model without bosonic 
degrees of freedom and with a dynamical symmetry breaking mechanism which needs a 
peculiar determination of the critical point. The main point of our simulation is that a not 
so critical tuning of its parameters gives a performance comparable to that of the Hybrid 
Monte Carlo. This is quite important with an eye on realistic models, like QCD, where 
a precise tuning is too expensive to be done. Moreover, the numerical precision argument 
makes us confident that the present algorithm can be valuable as an independent check 
for eventual biases of the Hybrid Monte Carlo. Simple scaling arguments (possibly naive 
because of the peculiar negation of momenta) indicate that the proposed algorithm should 
improve its performances with greater volumes or smaller fermion masses. 
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TABLES 



TABLE I. Horowitz e = 0.09 k = 6 



(5-0,7) 




TS2 


0.2) 


53(7) 


42(6) 


(-, 0.5) 


An/ r'\ 

49(5) 


40(4) 


(-, 1.0) 


44(4) 


37(3) 


(-, 2.0) 


o /i^\ 

58(7) 


46(5) 


(.-, 5-Oj 




4/(^0) 


(+, 0.2) 


43(5) 


35(4) 


(+, 0.5) 


38(3) 


29(3) 


(+, 1.0) 


44(5) 


31(3) 


(+, 2.0) 


41(4) 


29(2) 


(+, 5.0) 


66(9) 


41(5) 


TABLE II. Hybrid (rescaled by C) e = 0.06 N^d = 8 


So 


rs 






44(3) 


32(2) 


+ 


26(2) 


19(1) 
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